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SECTION  I 
INTRODUCTION 

During  the  last  decade,  many  successful  computational  methods  of  solution 
of  two-dimensional  and  axisynmetric  transonic  flows  have  appeared.  Most  early 
methods  (References  6 and  7)  depended  on  simple  geometry  and  the  small  pertur- 
bation potential  equation.  Jameson  (Reference  21)  was  one  of  the  first  to 
realize  the  utility  of  non-Cartesian  coordinate  systems  for  the  solution  of 
non-simple  geometry  with  the  exact  velocity  potential  equation.  The  present 
work  seeks  to  extend  this  idea  to  include  the  numerical  body-fitted  coordinate 
system  generation  method  developed  by  Thompson  (References  1 and  4).  Later  work 
by  Thames  and  Thompson  along  lines  similar  to  the  present  work  is  presented  in 
Reference  22. 

In  the  current  work,  a quick,  simple  method  of  generating  simply  connected 
coordinate  systems  is  developed  for  thin  and  slender  bodies,  while  Thompson's 
original  grid  generation  method  is  reserved  for  more  complex  bodies.  The  rrixed 
differencing  scheme  of  Murman  and  Krupp  (Reference  7)  is  applied  to  the  trans- 
formed exact  velocity  potential  equation.  The  method  is  verified  by  applying 
it  to  several  bodies  of  Interest. 
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SECTION  II 

NUMERICAL  GENERATION  OF  CONTROL  VOLUME  BOUNDARY  - 
FITTED  CURVILINEAR  COORDINATE  SYSTEMS 

For  the  numerical  solution  of  fluid  flow  problems,  it  would  be  convenient 
to  have  a curvilinear  coordinate  system  whose  outer  boundaries  are  coincident 
with  the  desired  solution  control  volume.  Thompson,  Thames,  and  Mastin  (Refer 
ence  1)  have  developed  such  a numerical  method.  Thames  (Reference  2)  provides 
perhaps  the  most  complete  and  comprehensive  explanation  of  the  method.  The 
first  part  of  this  discussion  is  based  on  the  work  of  Thames.  It  provides  a 
basic  explanation  of  the  method  and  the  changes  made  to  apply  it  efficiently 
to  the  solution  of  the  full  velocity  potential  equation  for  compressible  two- 
dimensional  and  axisymmetric  flows  in  a simply  connected  region.  A brief  expla- 
nation of  the  numerical  analysis  used  to  complete  the  transformation  is  presented. 

1.  BASIC  TRANSFORMATION 

Consider  transforming  an  arbitrary  two-dimensional,  simply  connected 

region,  D,  such  as  that  shown  in  Figure  1 onto  a rectangular  region,  D*.  as 

★ * * 

shown  in  Figure  2.  We  require  that  map  onto  r^,  onto  r^,  onto 
and  r4  onto  r^.  Region  D is  referred  to  as  the  physical  plane  and  region  D* 
as  the  transformed  plane.  For  the  present  purposes,  Tj  is  identified  as  a 
body  surface,  and  are  considered  to  be  straight  lines  of  symmetry  for 
a possible  axisymmetric  consideration!  is  taken  as  an  outer  boundary.  In 
the  present  work,  r3  will  be  used  as  both  an  "infinity"  semi-circle  and  a 
straight  wind-tunnel  wall.  Note  that  all  of  the  transformed  boundaries  are 
constant  coordinate  lines  ( n-lines)  in  the  transformed  plane. 

+ Note  that,  in  this  discussion,  no  generality  is  lost  by  this  imposed 
condition,  since  and  r2  are  treated  as  being  arbitrarily  shaped. 
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For  mathematical  reasons  discussed  by  Thames  (Reference  2)>Laplace 
equations  are  taken  as  the  generating  elliptic  system 


K + F =0 
xx  yy 


nxx  + nyy  = 0 


(1) 

(2) 


with  the  Dirichlet  boundary  conditions 


where  ni . £2»  n3»  and  F4  are  constant  and  ^(x.y),  n2^x »>" ) » ^(x.y).  and 
n4(x,y)  are  specified,  non-constant  functions  on  * r2>  r3,  and  r4» 
respectively. 


For  the  method  to  be  practical,  it  is  necessary  to  perform  all  numer- 
ical computations  in  the  uniform  rectangular  transformed  plane.  Following 
Thames  (Reference  2),  Equations  1 through  6 transform  to  the  system 


“XFF  - 2 6x?n  + yx^ 


ayFF  - 2eyCp  ♦ Yynn 


= 0 
= 0 


(7) 

(8) 


where 


2 2 
a=  K+y 
n Jn 


6=  xfx  + y,  y 
5 n JF  •'n 


2 2 
Y=XF  +y^z 


(9) 

(10) 

(ID 


with  the  transformed  boundary  conditions 

. [F  ,n  ] e r/ 
5 


■ "'W*,  WVUIIUUI 

R.KM 

□ Lf2(5>nl)J 


(12) 


= 

Fi(x.y) 
,ni  , 

, (x,y)  c Tj 

(3) 

w— 

h 

= 

F2 

n2(x,y) 

» a 

, (x,y)  e r2 

(4) 

F 

n 

■ a 

= 

F3 (x,y) 
n3 

, (x,y)  e r3 

(5) 

5 

n 

s 

nj(x.y) 

, (x,y)  e r4 

(6) 
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X 

9l 

y 

92(^2,n^ 

X 

h-|  U *93 ) 

y 

h2U» n3 ) 

X 

^ U4.n) 

y 

k2U4.11) 

» Cc.nJ  c r2 
» U.n]  g r3 
. U.n]  e r4* 


(13) 

(14) 

(15) 


The  functions  f ^ , f 2,  h-j,  and  h2  are  specified  by  the  known  shape  of  the 
boundaries  and  and  the  specified  distribution  of  £ thereon.  Simi- 
larly, g.| , g2>  k.| , and  k2  are  specified  by  the  desired  shape  of  the  boun- 
daries r2  and  r4  and  the  specified  distribution  of  n on  those  boundaries. 


The  set  of  Equations  7 through  15  is  a quasi -linear  elliptic  system 
for  the  physical  coordinates,  xU,n)  and  y(C,7i),  in  the  transformed  plane. 
With  this  system,  (x,y)  are  specified  for  each  U,n)  point  on  the  four 
boundaries.  Coordinates  U,n)  are  conveniently  numbered  as  1,  2,  3,  4, 
....  and  1,  2,  3,  . . .,  n^.  This  simple  numbering  is  the  con- 
stant step  size  in  the  computational  transformed  plane.  With  the 
boundaries  specified.  Equations  7 and  8 are  numerically  solved.  The 
resulting  natural  coordinate  system  has  an  n-line  coincident  with  the 
ri  and  boundaries  and  an  S-line  coincident  with  the  r2  and  r4  boundaries. 


Since  this  method  does  not  produce  orthogonal  coordinates,  the  fluid 
flow  equations  to  be  solved+  on  the  natural  coordinates  must  be  transformed 
by  direct  implicit  partial  differentiation.  Thames  (Reference  2)  gives 
the  derivative  expressions  for  time-dependent,  two-dimensional  systems, 
and  shows  that  the  type  (i.e.,  elliptic,  parabolic,  hyperbolic)  of  the 
fluid  equations  does  not  change  under  the  transformation. 

+Note  that  this  method  applies  to  any  set  of  partial  differential 
equations. 
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The  numerical  method  used  to  solve  Equations  7 and  8 was  based  on 
a modification  of  a computer  code  by  Thompson  (Reference  1)  to  solve 
only  simply  connected  regions  with  various  options  added  that  will  be 
discussed  later. 

As  discussed  by  Thames  (Reference  2),  a reasonable  "Initial  guess" 
(i.e.,  iterate  #0)  Is  required  to  guarantee  a converged  solution  for 
Equations  7 and  8 when  using  accelerated  successive  over  relaxation  (SOR) 
iteration.  Figures  3 and  4 are  samples  of  coordinate  system  Initial 
guesses  for  two  bodies  for  which  the  compressible  transonic  velocity 
potential  solutions  are  discussed  In  Section  IV.  Figure  3 is  a NACA  0012 
airfoil  with  an  outer  boundary  of  a circle  of  radius-ten  body  lengths, 
which  represents  an  "Infinity"  boundary.  The  r2  and  r4  boundaries  are 
held  straight  for  convenience  as  llnes-of-flow  symmetry.  The  interior 
of  the  field  is  a combination  of  straight  lines  and  ellipses.  Figure  4 
is  a numerical  model  of  a NASA  wind-tunnel  model  in  a simulated  NASA 
Langley  Research  Center  16T  wind  tunnel  used  in  an  AGARD  "Improved  Noz- 
zle Testing  Techniques"  study  (Reference  12). 

2.  EXTENSION  OF  THE  BASIC  TRANSFORMATION  TO  PROVIDE  FOR  CONTRACTION 
OF  THE  GRID 

The  adjustment  of  grid  spacing  may  be  necessary  in  areas  of  high  flow 
gradients.  Figures  5 and  6 demonstrate  that  the  solution  results  of  the 
elliptic  system  of  Equations  7 and  8 are  good;  however,  when  the  flow 
gradients  are  large,  the  grid  is  not  tight  enough  to  ensure  that  truncation 
error  does  not  affect  the  solution.  Also,  for  an  axi symmetric  problem, 
close  grid  spacing  near  the  axis  (y=0)  Is  necessary  for  an  accurate  solu- 
tion. More  control  can  be  obtained  over  both  5 and  n-Hne  spacing  by 
alteration  of  the  elliptic  system  of  Equations  7 and  8. 

Thompson,  et  al  (Reference  1)  have  shown  that  a set  of  equations  of  the 
form 


V2C  = p(x,y) 

(16) 

v2n  = Q(x,y) 

(17) 

can  be  formulated  with  the  boundary  conditions  given  In  Equations  3 thru  6. 
P = P(x,y)  and  Q = Q(x,y)  are  continuously  differentiable  negative 
functions  on  D and  its  boundary. 
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Study 


Figure  5.  NACA  0012  Airfoil  - Physical  Coordinate  System  with 
No  Contraction 
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Figure  6.  Agard  AD  HOC  Study  - NASA  Forebody  - 15°  Boattail 
NASA  Tunnel 

Physical  Coordinate  System  with  No  Contraction 
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By  interchanging  dependent  and  independent  variables  in  Equations 
16  and  17  we  can  derive  the  non-homogeneous  set 

aXEE  ‘ 2BxCn  + YXnn  = PU’n)  * xn  QU>n)]  (18) 

ay^  - 2By^  + yy^  = -J2[y^  P(C,n)  + yn  Q(C.n)]  (19) 

where  a,  6,  and  y are  as  defined  earlier  and  the  boundary  conditions  are 
given  by  Equations  12  through  15. 

Studies  by  Thompson,  et  al  (Reference  1)  have  demonstrated  that  the 
exponential  function,  because  it  is  harmonic,  is  a good  choice  for  the 
functions  P and  Q in  the  above  formulation. 

The  following  forms  for  P and  Q allow  attraction  to  specified  E 
and  n-lines  and/or  points. 


P(5,n)  = -Y  A,  Cjx,yj  ~ Ck  e 

l?(x’y)  ‘ ?kl 


-Ck  |E(x,y)  - Ek  ( 


m 

% 


B ^x,y)  " ^ t e 1 £ 
|E(x,y)  - Ep  ' 


(20) 


Qtt.n)  - -f  A'k  e 

|n(x,yj  - nk 


-Ck  ( n(x,y)  - nk 


y h 

k 


n(x,y)  - n 


l e 


»y)  - OjT 

where  R ="\/[E(x,y)  - E„]2  + [n(x,y)  - n„]2'  , 


(21) 


(22) 


and  n(n'  )is  the  number  of  attraction  E-lines  ( n-lines)  and  m(m  ) is  the 
number  of  attraction  ^-points  (n-points).  Note  that  E(x,y}/  Ek,E^  and 
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n(x,y)  } nk>  n^.  The  first  term  of  Equation  20  has  the  effect  of  attract- 
ing the  5-lines  to  the  5k-lines.  The  second  term  of  Equation  20  causes 
the  attraction  of  5-lines  to  the  points  [5^,  n^].  Similar  results  are 
obtained  with  the  two  parts  of  Equation  21  with  the  n-lines  to  the  n^-lines 
and  the  points  [5^,  n^]. 

The  attraction  concept  is  a useful  tool;  however,  its  application 
requires  trial  and  error  for  the  beginner  and  experience  in  order  to 
quickly  obtain  a suitable  coordinate  system. 

The  above  method  of  attraction  was  found  to  work  well  when  it  is 
desirous  to  attract  to  the  5k=l  or  the  nk=l  lines  only.  The  method  was 
also  suitable  for  contraction  to  a line  in  the  field.  However,  when  it 
was  found  to  be  necessary  to  attract  to  the  line  in  addition  to 

the  nk~l  line,  the  results  were  often  erratic.  This  attraction  to  both 
boundaries  was  found  to  be  necessary  to  prevent  the  grid  size  at  the 
outer  boundary  from  becoming  larger  than  1.0,  for  example,  in  the  case 
where  large  amounts  of  contraction  were  desired  at  the  body  (nk=l)  bound- 
ary. Instead  of  holding  the  grid  spacing  nearly  constant  at  the  outer 
boundary,  all  of  the  n-lines  tended  to  collapse  away  from  the  outer 
boundary  and  locate  in  the  bottom  part  of  the  field  near  the  area  of 
high  attraction. 

In  an  attempt  to  cure  the  problem,  the  forms  of  P and  Q have  been 
re-defined  as 


p(5,n)  = £ Ak  e 'Ckl?(x’y)  " Ck*+  £ 


Qfe.n)  - £ Ak  e 

k=l 


where  R£  is  still  defined  by  Equation  22.  This  modification  was  developed 
jointly  by  the  author  and  Hodge  (Equation  22).  The  coefficients  are  as 
defined  for  Equations  20  and  21.  For  attraction  to  a line  or  point  at  n=l 
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(or  £=1),  the  Ak's  and  B^'s  should  be  negative.  This  ensures  that  P and 

Q are  negative  in  the  area  of  the  lower  (i.e.,  £ or  n=l)  boundary.  For 

attraction  to  an  upper  boundary  (£=£  , or  n=n  , ) the  attraction  coeffi- 

max  max 

cients  should  be  positive. 

Equation  23  is  different  from  Equation  20  in  sign  and  the  deletion 
of  the  (£(x,y)  - £k)/|£(x,y)  - £|J  term  which  is  a SIGN  function.  It 
has  been  determined  that  this  SIGN  function  was  the  cause  of  the  erratic 
behavior  discussed  earlier.  Equation  24  was  changed  in  a similar  manner. 

Figures  7 and  8 are  plots  of  the  solution  of  Equations  18  and  19 
with  P and  Q defined  by  Equations  23  and  24.  The  coefficients  for  P and 
Q are  listed  in  Tables  1 and  2 for  Figures  7 and  8,  respectively.  The 
resulting  grids  are  obviously  more  suitable  to  numerical  calculation. 

The  grid  spacing  has  been  adjusted  to  account  for  the  flow  gradients  in 
the  velocity  potential  flow  solution  described  in  Section  III.  The 
boundary  conditions  can  be  applied  on  coordinate  lines  which  are  coinci- 
dent with  the  boundaries.  These  two  coordinate  systems  have  been  used 
successfully  to  solve  for  the  velocity  potential  for  flows  which  are 
totally  subsonic. 

3.  OPTIMIZING  THE  COORDINATE  GENERATION 

In  general,  the  convergence  rate  for  the  numerical  solution  of  a 
partial  differential  equation  solved  on  an  orthogonal  grid  is  faster 
than  one  solved  on  a non-orthogonal  grid.  With  an  orthogonal  grid,  the 
coefficients  of  the  off-diagonal  terms  are  minimized.  With  this  in  mind, 
several  attempts  have  been  made  to  make  the  grid  system  more  orthogonal 
by  altering  the  equation  set  18,  19,  and  12  through  15. 


The  first  procedure  merely  involved  changing  the  boundary  conditions 
* * 

on  and  from  Dirichlet  to  Neumann  expressions.  The  requirement  for 
orthogonality  from  Equations  18  and  19  is 


Figure  7.  NACA  0012  Airfoil  - Physical  Coordinate  System  with 
Contraction 
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Along  and  r4,  y =0  for  the  physical  coordinate  systems  illustrated 
thus  far.  Therefore,  Equation  25  reduces  to 

Vn  = 0 

However,  as  x ^0,  we  must  have  x =0.  Therefore,  the  new  boundary  condi- 
* n * £ 

tions  for  r2  and  f4  are  as  follows: 


V 

‘o 

y 

~^2  ^2 

V 

'o 

y 

^(£4  »h) 

[?,n]  e 

C5.n]  e rj 


(26) 

(27) 


Figure  9 is  a solution  of  the  equation  set  18,  19,  12,  26,  14,  and  27. 

★ * 

For  the  first  iterate,  x must  still  be  specified  along  and  r^.  For 

all  iterations  after  the  first  one,  the  value  of  x is  calculated  from 

a one-sided  second-order  difference  expression.  As  explained  by  Thames 

(Reference  2),  the  use  of  Dirichlet  boundary  conditions  with  Equations 

16  and  17  guarantees  that  the  solution  satisfies  the  maximum  principle,  which, 

in  turn,  guarantees  that  all  extrema  occur  on  the  boundaries  of  region  D. 

However,  with  the  use  of  the  Neumann  boundary  condition  shown  In  Equations 

26  and  27,  the  above  is  no  longer  the  case.  The  problem  which  can  result 

is  illustrated  in  Figure  10,  which  is  an  enlargement  of  the  rear  portion 

of  the  body  of  Figure  9.  The  points  (1,2)  and  (?^v*2)  both  are  outside 

the  boundaries  specified  by  F^.  This  problem  was  solved  by  applying  the 

boundary  conditions  of  Equations  26  and  27  at  each  Individual  point  on 
★ ★ 

r2  and  only  as  long  as  the  x-distance  between  points  did  not  become 
less  than  a specified  input  value.  The  results  of  such  a requirement  are 
illustrated  in  Figures  11  and  12,  where  the  smallest  step  size  allowed  was 
0.005  and  0.00375,  respectively.  Obviously,  the  Neumann  boundary  condi- 
tion is  not  valid  at  the  points  where  the  minimum  step  size  was  enforced, 
but  the  results  are  still  excellent. 
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9.  Agard  AD  HOC  Study  - NASA  Forebody  - 15°  Boattail 
NASA  Tunnel 

Physical  Coordinate  System  with  Contraction  and 
Neumann  Boundary  System 


Figure  11.  NACA  0012  Airfoil  - Physical  Coordinate  System  with 
Contraction  and  Neumann  Boundary  Condition 


Agard  AD  HOC  Study  - NASA  Forebody  - 15  Boattail 
NASA  Tunnel 

Physical  Coordinate  System  with  Contraction  and 
Controlled  Neumann  Boundary  Condition 
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Table  3 is  a list  of  the  convergence  requirements,  field  size,  number 
of  Iterations  required  for  convergence,  and  CDC  6600  computation  time 
required  for  convergence  for  all  of  the  coordinate  systems  of  Figures  7 
through  12  and  Figures  15  and  17,  described  in  the  following  sections. 

The  data  list  includes  the  time  required  to  create  the  initial  specifi- 
cation, the  solution  of  Equations  7 through  15,  and  the  solution  of 
Equations  18  and  19  with  the  appropriate  boundary  conditions  (i.e., 
Dirichlet  or  Neumann).  The  number  of  iterations  required  for  convergence 
is  only  that  number  used  in  the  satisfaction  of  Equations  18  and  19. 

As  can  be  seen  from  Table  3,  the  time  required  to  compute  the  very 

non-orthogonal  system  of  Figures  7 and  8 is  a barely  acceptable  (i.e., 

for  the  compressible  velocity  potential  solution)  31  to  41  seconds.  To 

★ ★ 

apply  the  orthogonal  boundary  condition  along  ^ and  involves  more 
than  doubling  the  computation  time.  Clearly,  a more  efficient  method  of 
producing  acceptable  coordinates  is  needed  for  use  with  a velocity  poten- 
tial solution. 


TABLE  3 

COORDINATE  SYSTEM  COMPUTATIONAL  SUMMARY 


Figure 

Number 

Convergence 

Requirement* 

Field  Size 
^max  x nmax^ 

Number  of 
Iterations 

Time 

(Seconds) 

7 

37  x 30 

80 

31 

8 

0.0001 

61  x 25 

72 

41 

9 

0.0001 

61  x 25 

259 

83 

11 

0.0001 

37  x 30 

351 

73 

12 

0.0001 

61  x 25 

257 

84 

15 

0.00027 

81  x 15 

350 

75 

17 

0.0001 

111  x 20 

65 

23 

*For  all  points  in  the  field,  the  largest  change  from  one  iteration  to 
the  next  in  both  x and  y is  less  than  the  specified  value. 
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4.  DEVELOPMENT  OF  A STREAMLINE  TYPE  COORDINATE  SYSTEM 

The  coordinate  systems  developed  in  Subsection  C are  very  useful  for 
the  solution  of  the  velocity  potential  equation  as  long  as  the  flow  is 
totally  subsonic  (i.e.,  the  equation  type  is  elliptical).  As  the  flow 
becomes  supersonic  the  equation  type  is  hyperbolic  and  one-side  upwind 
differences  must  be  employed  for  computational  stability.  The  problem 
with  the  coordinate  systems  discussed  until  now  is  that  "upwind"  is  not 
readily  evident.  This  problem  will  be  discussed  fully  in  Section  III. 

One  solution  to  the  problem  is  the  introduction  of  a streamline  type 
of  coordinate  system.  Figure  13  illustrates  the  new  system.  The  system 
is  mainly  oriented  toward  the  "wind-tunnel"  type  of  boundaries.  With  this 
system  the  n=constant  lines  are  obviously  in  the  streamline  direction  and 
allow  the  easy  application  of  one-sided  upwind  derivatives  for  the  flow 
solution. 


Figure  14  is  the  initial  specification  of  the  ogive-cylinder-boat- 
tail  body  of  Figure  12.  In  Figure  14,  a new  initial  specification  method 
was  used  which  was  suggested  by  Ghia  (Reference  3).  The  method  has  not 
yet  been  adapted  to  the  coordinate  system  of  Figure  1.  (For  future  refer- 
ences, the  coordinate  system  of  Figure  1 will  be  referred  to  as  the 
"elliptical"  type  of  coordinate  system  and  that  of  Figure  14  as  the 
"streamline"  type).  Figure  15  represents  the  best  coordinate  system 
found  with  this  procedure.  Figure  16  Is  a plot  of  P(C,n)  used  in  the 
calculation.  Neumann  boundary  conditions  were  applied  on  r^,  r3>  and 
r4.  This  coordinate  system  was  used  to  solve  problems  which  had  super- 
sonic "pockets"  in  the  boattail  area  with  good  success.  However,  as  can 
be  seen  from  Table  3,  the  72.7  sec  computation  time  for  the  creation  of 
the  coordinate  system  was  unacceptable. 
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Figure  14.  Agard  AD  HOC  Study  - NASA  Forebody  - 15  Boattail  - 
NASA  Tunnel 

"Streamline"  Type  Physical  Coordinates  System  Initial 
Specification 


Agard  AD  HOC  Study  - NASA  Forebody  - 15  Boattail 
NASA  Tunnel 

"Streamline"  Type  Physical  Coordinate  System  with 
Contraction 
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To  solve  this  computation  time  problem  and  to  achieve  more  orthogonal 
systems,  the  following  method  suggested  by  Ghia  (Reference  3)  was  imple- 
mented. Consider  writing  Equations  18  and  19  as 


(xntyn>  x«  ' 2(xeV  Vn)x5n  * (xH>  x„n 


+ J2(x^P  + xnQ)  = 0 


(xn*,n)  - 2IW  y5Vy5n  * %, 


+J2(ycP  + ynQ)  = o 


Consider  a boundary  such  as  r2  and  r4  of  Figure  13  which  are  straight  and 
aligned  with  the  y-^£.  In  this  case,  x^O.  If  one  desires  an  orthogonal 
system,  y^=0  (1  .e. , the  line  segment  for  at  least  the  first  three  points, 
£=1,  2,  and  3,  is  straight)  and  y^=0  Now  applying  these  to  Equation  29, 

we  find  x2ynn  + J2ynQ  = 0,  but  J = x?yn  - xpy?  = x^.  Therefore, 

q - -y„n/y2  (30) 


Which  states  that  the  value  of  Q necessary  to  have  an  orthogonal  system 
near  the  boundaries  r2  and  can  be  simply  calculated  from  the  specified 
distribution  of  points  along  that  boundary.  To  evaluate  Q over  the 
interior  of  the  field,  make  a linear  distribution  of  Q between  the  two 
boundaries  for  each  rrline,  i.e.. 


Q(C,n)  Q(£|^x>n)  - ,n)j 


CMAX  " C 


30 


OTMMMi 
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Equation  30  has  been  found  to  be  easily  applied  and  very  effective. 
If  a third  order  polynomial  is  used  to  specify  the  initial  distribution 
of  y along  and  r^,  the  four  coefficients  of  the  polynomial  can  be 
determined  by  specifying  y,  y^,  and  y at  n=1  (the  lower  boundary)  and 

y at  n=nma„  (the  upper  boundary).  Normally,  y is  chosen  to  be  zero 

maX  nn 

and  yn  is  the  desired  step  size  at  the  lower  boundary.  The  n-line 
distribution  of  Figure  17  is  governed  by  Q determined  from  Equation  30. 
The  input  value  of  y^  at  the  lower  boundary  was  0.0095.  The  converged 
solution  had  an  average  step  size  of  0.0097  between  the  n=l  and  n=2 
lines. 


An  analogous  equation  can  be  derived  from  Equation  28  for  the 
boundary  r3  of  Figure  13  when  it  is  straight  and  aligned  with  the  x-axis, 
that  is. 


(32) 


At  the  upper  boundary  of  Figure  17,  P was  determined  from  Equation  32. 


An  expression  for  P can  be  derived  from  Equation  28  for  an  arbitrary 

boundary  such  as  of  Figure  13.  Assume  that  along  this  coordinate 

line  (n=l)  xnn=0-  This  means  that  the  £-lines  are  straight  for  the  small 

segment  n=l»2,3.  As  the  end  boundaries  and  r«  are  specified,  y 

may  be  calculated  at  n=l  at  each  end  (£=1  and  €s£maw).  The  value  is 

max 

calculated  by  a second-order  one-sided  difference  expression,  since 
this  Is  the  value  used  in  the  numerical  calculations  and  not  the  input 
value  used  in  the  polynomial  discussed  above.  A linear  distribution 
of  y between  y (1,1)  and  y (£m,„,l)  is  taken  and  held  constant  for 
the  calculation  of  P.  Similarly,  y^  may  be  determined  along  r-j  from 
the  specified  values  of  y at  the  end  points.  Then  x , y , x^  , and 
y may  be  calculated  from  central  difference  expressions  from  the  input 
specifications  of  x and  y along  r^.  Since  an  orthogonal  system  is 
desired,  and  B=0  is  the  condition  for  orthogonality,  x^  may  be  determined 
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Figure  17.  Agard  AD  HOC  Study  - NASA  Forebody  - 15°  Boattall 
NASA  Tunnel 

"Streamline"  Type  Physical  Coordinate  System  with 
Calculated  Contraction  Factor 
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along  r1  from  6=0,  i.e.,  x^  -y^/x^-  The  Jacobian,  J=x^yn  -x^,  is 
now  known.  Q is  determined  from  Equation  30  at  S=1  and  and  a 

max 

linear  distribution  between  the  two  values  is  taken.  From  the  above 
assumption  and  now  known  values  along  1^,  the  following  equation  can  be 
applied: 


The  very  successful  results  are  shown  in  Figure  17.  Because  Q usually 
varies  very  rapidly  in  the  vicinity  of  n=l,  experience  has  shown  that 
the  Q used  in  Equation  33  should  be  an  average  of  the  values  at  n=l 
and  n=2. 


Thompson  (Reference  4)  showed  that  a truncation  error  resulting 
in  a numerical  diffusive  effect  occurs  in  a flow  calculation  when  the 
second  derivatives  of  the  physical  coordinates  (x^  in  this  case)  are 
large  in  regions  where  the  dependent  variables  have  significant  second 
derivatives  in  the  direction  normal  to  closely  spaced  coordinate  lines. 
The  use  of  the  above  method  for  determining  P brought  this  problem  to 
light.  The  method  forces  the  grid  spacing  specified  at  the  wall  to  be 
continued  through  the  field.  For  example,  large  values  of  ^ at  the 
wall  forced  the  problem  to  propagate  throughout  the  field.  The  result 
in  the  potential  flow  calculations  was  that  the  flow  would  not  accelerate 
to  the  expected  values  in  regions  of  large  gradients  in  the  flow  direc- 
tion. It  was  determined  by  trial  and  error  and  comparison  with  experi- 
mental and  theoretical  data  that  the  condition 


|x^|  < 0.2|x?| 

was  sufficient  to  ensure  minimum  influence  on  the  flow. 


(34) 


33 
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SECTION  III 

APPLICATION  TO  THE  EXACT  VELOCITY  POTENTIAL  EQUATION 

The  natural  coordinate  system  method  developed  in  Section  II  is 
applied  to  the  numerical  solution  of  the  exact  compressible  velocity 
potential  equation  for  two-dimensional  and  axisymmetric  bodies  in  a 
simply  connected  field.  The  velocity  potential  equation  is  presented 
in  cartesian  form  and  then  as  it  appears  in  the  transformed  plane.  The 
boundary  conditions  are  presented  in  both  forms.  The  numerical  method 
of  solution  is  explained.  The  method  of  integrating  the  pressure  distri- 
bution to  obtain  boattail  drag  on  axisymmetric  bodies  is  developed. 

1.  THE  COMPRESSIBLE  VELOCITY  POTENTIAL  EQUATION 

For  steady,  invicid,  irrotational , isentropic,  compressible  flow, 
the  exact  velocity  potential  equation  in  non-dimensional  form  is 

Ad>  + B<f>  + C<p  + Dd>  = 0 
^xx  Txy  yyy  y 

where  the  coefficients  are  defined  as 


A ■ - »x 

(36a) 

B = -2d>  <b 

Vy 

(36b) 

2 2 

C - a - »y 

(36c) 

D = cta2/y 

(36d) 

Two-dimensional  flow  is  represented  when  a=0  and  axisymmetric  flow 
when  a=l.  The  local  speed  of  sound,  a,  is  defined  as 
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A variable  of  major  interest  in  boattail  performance  calculations  is  the 
pressure  coefficient  and  is  defined  from  Liepman  and  Roshko  (Reference  5) 
as 


C 


P 


_\  2+(y  -1)M2  , 

a 


(39) 


The  transformation  of  the  above  equations  into  the  computational 
(Cfn)-plane  is  a straightforward  process  with  the  use  of  Thames'  (Refer- 


ence  2)  thesis.  Equation  35  becomes 

kl%  + k2<J)Cn  + k3^nn  + k4*5  + k5<f>n  = 0 

(40) 

where 

lr»v"Wc,5 

(41a) 

k2=-[2Aycyn-B(x5yn*xnyEK2x5xn] 

(41b) 

k3=Ay|-Bxeye+cx2 

(41c) 

k4=-JCA(fliXn.Ajyn)tBB,*c(c,xn-c2yn)-J2°xn] 

(41d) 

k5=_3CA(A2yrAl  *E  )+BB2+C(C2yc'Cl  x£  )+j2°x^  J 

( 41  e ) 

and 

AryA 

(42a) 

2 2 

A,=y  xrr-2yry  xr  +ytrx 

2 'n  EE  E n En  E on 

(42b) 

VJKyr\V*(VnJ5-xcVn> 

(42c) 

l2,JlV!S-Vcn)*(YtVxa,(JE) 

(42d) 

2 2 

C,  = v y -2xrx  y.  +xt:  y 

1 n EE  E n En  E nn 

(42e) 

2 2 

C,=xr,  xrr-2xrx  xr  +xrxnn 

2 n EE  E n En  E nn 

(42f) 
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Equation  37  transforms  to 


V(ynVVn)/J 


V(xeVxA)/J 


(43a) 

(43b) 


2.  BOUNDARY  CONDITIONS 


Various  boundary  conditions  are  presented  below  which  can  be  combined 
to  solve  a problem  of  choice. 


a.  Uniform-Stream  Boundary 

u=u» 

v=0 


(44a) 

(44b) 


This  condition  is  easily  expressed  as 


»=x 


if  x +y 


(45) 

(46) 


b.  Solid  Wall  Boundary 
u=u(x,y) 
v=0 


(47a) 

(47b) 


By  utilizing  Equation  43b,  the  transformed  expression  for  this 
boundary  condition  Is 

We/Xe  (48) 


c.  Tangential  Flow  Boundary 

v-n  = 0 


(49) 


where 


v= 


V+V 


and 


(50) 
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Equation  50  is  the  expression  for  the  unit  normal  to  a line  of  constant 
n given  by  Thames  (Reference  2).  The  performance  of  the  dot  product 
indicated  by  Equation  49  results  in  the  transformed  expression 

V8VY  (51) 

Note  that  Equation  51  reduces  to  Equation  48  along  a constant 
n-line  which  Is  straight  and  parallel  to  the  x-axIs  (i.e.,  a boundary 
where  y^=0). 

d.  Stagnation  Point 

Generally,  it  would  be  desirous  to  impose  a complete  stagnation 
condition  (i.e.,  u=0,v=0)  at  the  leading  and  trailing  edges  of  bodies 
such  as  the  NACA  0012  airfoil  discussed  in  Section  II.  However,  an 
attempted  application  of  both  conditions  when  dealing  with  the  potential 
equation  results  in  an  over  specification  of  the  boundary  condition. 

Therefore,  the  boundary  condition  that  was  used  is 

v=u*i  (52) 

It  was  found  to  be  necessary  to  specify  u*  initially  as  u^  when 
the  flow  was  started  from  a uniform  stream.  The  values  of  u*  were 
gradually  reduced  to  zero  as  the  solution  was  iterated.  The  transfor- 
mation of  Equation  52  results  in  the  expression 

V(ynVJu*)/y5  (53) 

3.  NUMERICAL  METHOD 

In  1970,  Murman  and  Cole  (Reference  6)  and  Murman  and  Krupp  (Reference 
7)  demonstrated  that  it  is  necessary  to  use  a mixed  finite  difference 
scheme  to  solve  the  potential  equation  when  the  flowfleld  has  regions 
of  subsonic  and  supersonic  flow.  In  the  subsonic  (elliptic)  regions. 
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centra]  difference  expressions  are  necessary  to  ensure  numerical  stability. 
However,  in  supersonic  (hyperbolic)  regions,  one-sided  upwind  differences 
are  necessary. 

The  implementation  of  this  method  in  this  report  has  been  simple 
and  straightforward.  First,  the  appropriate  boundary  condition  is 
solved  along  the  constant  5=1  line  (r2  boundary).  For  the  "elliptical" 
type  of  coordinate  system  of  Figure  1,  the  tangential  flow  boundary 
condition.  Equation  51,  is  solved.  Both  an  explicit  method  and  an 
implicit  tri-diagonal  solver  have  been  used  successfully.  For  the 
"streamline"  type  of  system  in  Figure  14,  the  uniform  stream  boundary 
condition.  Equation  45,  was  implemented.  The  solution  then  marched  to 
the  C=2  line.  For  each  point  on  the  line  the  type  of  Equation  39  is 
determined  from  the  characteristic  equation  of  Equation  35 

A*x  - B*x  * C ■ 0 • (54) 

Specifically,  if 

B2-4AC<0  it  is  elliptic 

B2_4AC=0  it  is  parabolic 

B2-4AC>0  it  is  hyperbolic 

Thames  (Reference  2)  demonstrated  that  this  type  of  the  differential 
equation  is  invariant  under  the  coordinate  transformation.  A,B,  and  C 
are  determined  by  using  only  central  differences.  When  Equation  35  is 
either  elliptic  or  parabolic,  second  order,  central  derivatives  are 
used  to  determine  the  coefficients.  Equations  40  and  42  and  to  express 
Equation  39.  When  the  equation  type  is  hyperbolic,  first  order,  one- 
sided upwind  (decreasing  5)  differences  in  the  ^-direction  and  second 
order,  central  differences  in  the  n-di recti  on  are  used  to  determine 
the  coefficients  and  express  Equation  39.  The  appropriate  boundary 
conditions  are  expressed  in  a similar  manner  and  the  entire  £-line  is 
posed  in  tri-diagonal  form.  This  method  is  essentially  line  SOR. 
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After  each  line  is  solved  the  solution  is  weighted  with  acceleration 
parameters  which  have  different  values  for  the  elliptic  region,  the 
hyperbolic  region,  and  the  boundaries.  The  values  of  these  parameters 
are  discussed  in  Section  IV. 

Each  £-line  is  solved  in  turn  until  the  rear  boundary,  r4,  is 
reached.  It  is  solved  in  the  same  fashion  as  the  boundary.  This 
method  is  then  iterated  over  the  field  until  the  maximum  change  of  the 
potential  from  one  iteration  to  the  next  reaches  some  convergent  criteria. 

This  method  has  worked  well  for  subsonic  free-stream  flow  having 
embedded  supersonic  pockets  with  mild  shocks.  Thus  far,  the  maximum 
free-stream  Mach  number  that  has  been  successfully  run  is  0.99.  The 
"elliptical"  type  of  coordinate  system  does  not  work  when  a supersonic 
area  falls  in  a portion  of  the  grid  where  the  n-lines  are  not  well 
aligned  with  the  flow. 

4.  AXISYMMETRIC  PRESSURE  DRAG  INTEGRATION 

An  important  result  of  an  axisymmetric  potential  flow  solution  is 
the  integrated  pressure  drag.  From  Figure  18,  we  can  see  that 

P=  -Pn(n)  (55) 

where  nvl'  is  the  unit  normal  defined  outward  toward  increasing  n by 
Equation  50.  With  an  increment  in  arc  length,  ds  = /fdx)2+(dy)2,  giving 
an  increment  in  surface  area,  dS=ZIIyds,  the  drag  increment  is  defined  as 

dD=(P  • i )dS  (56) 

where  i is  a unit  normal  in  the  positive  x-di recti  on.  By  carrying  out 
the  operation  indicated  by  Equation  56, 

dD  = Py^dS/ /T  (57) 
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Note  that  the  normal  inclusion  of  the  step  size,  h,  in  the  integration 
formula.  Equation  65,  is  not  necessary  as  h=A£=An=l  by  definition  when 
the  coordinate  transformation  is  applied.  Such  is  the  case  in  all 
difference  formulas  used  in  this  report.  To  start  and  end  the  integration, 
the  trapezoidal  rule  is  used 
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SECTION  IV 
CALCULATIONS 

The  velocity  potential  solution  method  discussed  in  Section  III  was 
programmed  in  Fortran  IV.  In  this  section  the  basic  capabilities  and 
limitations  of  the  method  are  verified.  The  efficiency  of  the  numerical 
scheme  and  the  effects  of  acceleration  parameters  are  discussed.  The 
pressure  drag  integration  method  is  demonstrated.  Finally,  the  effects 
of  various  boundary  conditions  on  the  flow  over  a cone-cyl inder-boattail 
are  examined. 

1.  VERIFICATION 

A sphere  in  an  incompressible  flow,  a NACA  0012  airfoil  at  zero 
angle  of  attack,  an  ogive-cyl inder-boattail , and  a parabolic  arc  body 
of  revolution  were  chosen  for  basic  verification  of  the  method. 

Figure  19  shows  pressure  coefficient  contours  and  local  total 
velocity  contours  for  the  sphere.  From  Karamcheti  (Reference  8),  the 
minimum  pressure  coefficient  on  the  body  should  be  Cp=-1.25.  The 
calculated  value  of  Cp=  -1.244  is  0.48%  in  error.  The  calculated 
maximum  velocity  on  the  body,  Q=1.498,  is  very  close  to  the  ideal  value 
of  1.500.  The  "elliptical"  coordinate  system  used  to  calculate  the 
sphere  flow  was  essentially  a polar  one  with  the  far  field  at  15  body 
diameters  away. 

Figure  20  Illustrates  the  calculation  for  flow  over  a NACA  0012 
airfoil  at  zero  angle  of  attack.  The  data  comparison  for  the  incom- 
pressible case  Is  from  Abbott  and  Von  Doenhoff  (Reference  9).  The 
M=0.63  and  M=0.725  cases  compare  very  well  with  the  wind  tunnel  results 
of  References  9 and  10.  The  coordinate  system  used  was  that  of  Figure 
11  which  has  an  "infinity"  boundary. 
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To  check  supercritical  flowfield  calculations,  the  body  of  Figure 
17  was  run  at  Mach  numbers  of  0.4,  0.6,  0.8,  and  0.9.  The  results  of 
Figure  21  show  comparisons  with  wind-tunnel  data  from  References  12 
and  14  and  the  calculations  of  the  Douglass-Neumann  program  (Reference 
13).  The  results  are  quite  good  for  all  Mach  numbers.  The  solution  in 
the  shock  area  in  the  M=0.9  case  is  very  close  to  the  experimental  data. 

This  is  typical  of  non-conservative  formulation  of  the  governing  equations. 
Figure  22  is  a plot  of  the  pressure  distribution  for  the  whole  body. 

The  symbols  demonstrate  the  variable  grid  spacing  typically  used. 

The  coordinate  system  of  Figure  17  represents  the  NASA  LaRC  16T 
wind  tunnel.  The  wall  height  was  proportioned  to  1.5  body  lengths  as 
reported  by  NASA  LaRC  (Reference  12).  The  wall  boundary  condition  was 
simulated  with  a solid  wall.  Equation  48.  The  body  was  also  run  with 
the  upper  boundary  placed  at  5.0  body  lengths  away.  With  this  configuration, 
the  uniform-stream  boundary  condition.  Equation  45,  was  implemented. 

For  the  range  of  Mach  numbers  considered,  no  appreciable  effects  of  the 
upper  boundary  condition  were  seen.  Because  of  the  results  of  these 
calculations,  it  is  felt  that  the  solid-wall  boundary  condition  is 
sufficient  to  represent  the  real  LaRC  16T  slotted  wall  for  axisymmetric 
flows  with  only  small  supersonic  pockets  and  when  the  body  length  to 
test  section  radius  ratio  is  1.0  or  less. 

To  Investigate  flows  near  Mach  1.0,  the  10%  thick  parabolic  arc 
body  of  revolution  tested  by  Taylor  (Reference  15)  was  considered.  The 
80-inch  model  was  tested  in  the  NASA  Ames  14-Foot  Transonic  Wind  Tunnel. 

The  wind  tunnel  has  a square  ventilated  test  section  consisting  of 
longitudinal  slots  with  corrugated  inserts.  A turbulent  boundary  layer 
was  induced  on  the  model  at  R„  =27  X 106.  The  data  were  not  corrected 
for  wall  interference  effects.  The  body  was  truncated  at  X/L=0.854  to 
permit  mounting  on  the  sting. 
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Figure  21.  Pressure  Distributions  for  Flow  Over  the  Agard 

AD  HOC  Study  15°  Boattail  for  M=0.4,  0.6,  0.8,  and  0.9 


Figure  21.  (Concluded) 
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A calculation  was  made  for  this  body  at  M=0.975.  A solid  wind-tunnel 
wall  was  placed  at  1.0125  (81  inches)  body  lengths  radius.  The  results 
of  the  calculation  are  shown  In  Figure  23.  Also  shown  are  the  calculations 
of  Bailey  (Reference  16).  The  physical  step  size  along  the  body  was  a 
constant  0.015  for  0.104<.X/L<0.854.  Near  the  body  nose  the  step  size 
was  reduced  to  0.0065.  As  can  be  seen  from  Figure  23,  the  calculated 
solution  over-shoots  the  experimental  data  and  Bailey  results  over  the 
first  15%  of  the  body  as  if  the  nose  were  slightly  blunted.  Initially, 
this  problem  was  thought  to  be  caused  by  the  selection  of  x-step  size 
in  nose  area.  Numerous  variations  of  xU),  x'U)  and  x"U)  were  tried 
in  the  nose  area  to  study  the  grid  selection  effect  on  the  solution 
in  this  area.  The  results  were  always  the  same.  Refinement  and 
smoothing  of  the  grid  merely  resulted  in  better  definition  of  the  same 
solution. 

Another  possible  source  of  the  problem  is  the  manner  in  which 
coordinate  system  ^-derivatives  (i.e.,  x^,  x^,  y^,  y^,)  are  evaluated 
at  the  nose  on  the  lower  boundary.  For  the  solution  of  Figure  23, 
second-order  central  derivatives  are  used.  It  would  seem  apparent 
that  this  procedure  would  have  the  effect  of  smoothing  the  discontinuity 
at  the  nose  into  a cusp  with  the  expected  effect  of  causing  an  under 
expansion  in  the  nose  area.  However,  as  Figure  23  illustrates,  this 
was  not  the  case.  Note  that  as  the  tangency  condition.  Equation  49, 
is  solved  along  the  lower  boundary,  r-j,  and  not  the  potential  equation. 
Equation  39,  the  actual  value  of  y at  the  nose  point  does  not  enter  into 
the  calculations.  The  solution  along  r-j  was  also  separately  attempted  with 
second-order  upwind  and  downwind  differences  to  evaluate  the  derivatives 
at  the  nose.  In  both  cases,  the  solution  diverged  at  the  nose. 

The  calculated  solution  over  the  remainder  of  the  body  is  in 
excellent  agreement  with  the  experimental  data  of  Taylor  (Reference  15). 

The  solution  converges  well  In  all  areas  except  in  the  shock  region. 

In  the  150  iterations  allowed  In  the  calculation,  the  maximum  change  in 
$ at  the  last  Iteration  was  0.000012.  Elsewhere  in  the  field,  the 
convergence  was  better  than  0.0000005.  The  inflection  in  the  pressure 
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Figure  23.  Pressure  Distribution  and  Sonic  Line  for  M=0.975  Flow 
over  a 10%  Thick  Parabolic  Arc  Body  of  Revolution 
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distribution  curve  at  X/L=0.6  is  caused  by  the  presence  of  the  cylindrical 
sting.  Also  shown  in  Figure  23  is  the  sonic  line  in  the  flowfield.  The 
inflection  of  the  Cp  curve  is  also  present  in  the  sonic  line  as  the 
disturbance  propogates  almost  vertically  in  the  flowfield.  Calculations 
made  without  the  sting  did  not  show  an  inflection  in  either  the  surface 
pressure  distribution  or  the  sonic  line. 

2.  COMPUTATIONAL  EFFICIENCY 

The  calculations  of  Section  A provide  a data  base  to  examine  the 
efficiency  of  the  method.  When  programmed  on  the  CDC  6600  FTN  (0PT=1 ) 
compiler,  the  calculation  rate  for  transonic  flows  is  3700  grid  points 
per  second.  For  totally  subsonic  flows,  where  the  extra  supersonic 
coding  is  deleted,  the  computation  rate  increases  to  approximately  5000 
grid  points  per  second.  The  core  storage  for  a 101  x 22  field  size, 
such  as  that  used  for  the  results  of  Figure  23,  is  154000B.  The 
computation  time  for  axi symmetric  flowfields  of  this  type  is  30  seconds 
to  3 minutes,  depending  on  the  starting  conditions  and  the  closeness 
to  a free-stream  Mach  number  of  unity.  Two-dimensional  calculations 
characteristically  require  more  iterations  with  1 to  5 minutes  of 
computations. 

The  effects  of  acceleration  parameters  on  the  convergence  of 
various  flow  problems  is  illustrated  in  Figure  24.  Case  I is  the  only 
two-dimensional  condition  considered  due  to  the  larger  costs  of  2-D 
solutions.  The  primary  investigation  was  made  holding  ABC=0.85  while 
AER  was  varied  from  a low  of  1.1  to  a high  of  1.5.  When  AER  was  raised 
above  1.5  the  solution  diverged,  no  matter  what  the  value  assigned  to 
ABC.  To  study  the  effect  of  ABC  at  optimum  AER,  ABC  was  varied  from 
0.65  to  l.tf.  With  ABC  above  1.65,  the  calculations  diverged.  As  can 
be  seen  from  Figure  24,  with  the  optimum  acceleration  parameters  of 
AER=1 . 5 and  ABC=1.65,  NIC  could  be  reduced  to  355  with  103.1  seconds 
of  CDC  6600  computation  time.  The  worst  case  shown,  AER=1.1,  ABC=0.85, 
required  544  iterations  and  157.6  seconds  for  an  increase  in  time  of 
53%.  The  selection  of  optimum  acceleration  parameters  is  well  worth 
the  effort. 
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Case  II  is  an  almost  completely  subsonic  axisymmetric  cone-cyl inder- 
boattail  configuration.  AER  was  varied  from  0.85  to  1.8  and  ABC  was 
varied  from  0.6  to  1.1.  Changing  AER  causes  the  expected  trend  in  NIC; 
however,  varying  ABC  did  not  affect  the  solution  at  all  until  it  suddenly 
diverged  when  ABC  was  increased  to  1.0  from  0.95. 

Cases  III,  IV,  and  V are  different  conditions  solved  on  the  same 
body.  In  all  three  cases,  ABC  had  no  effect  on  the  convergence  until 
divergence  occurred  at  a value  of  ABC  near  unity.  Another  problem 
which  must  be  considered  is  that  at  the  lower  values  of  AER  the  solution 
will  not  develop  as  much  as  it  will  when  accelerated  near  the  limit. 

The  dashed  portions  of  the  lines  represent  these  areas.  The  less  dis- 
turbed the  field  is  (i.e. , over  a smooth  axisymmetric  body),  the  greater 
the  problem.  In  Case  I,  the  maximum  discrepency  in  Cp  on  the  body  sur- 
face was  only  ACp=0.001;  but  for  Cases  II,  III,  and  IV,  the  value  reached 
0.009.  To  ensure  the  most  correctly  developed  solutions,  it  is  necessary 
to  use  a value  of  AER  within  0.2  of  the  maximum  allowable. 

Case  V presents  an  especially  troublesome  case,  in  which  NIC 
increases  as  AER  is  increased.  This  phenomenon  can  be  credited  to  the 
lower  AER  values  damping  out  the  small  change  in  the  flowfield  when  an 
M=0.9  solution  is  started  from  an  M=0.8  solution.  The  calculation 
savings  in  obtaining  the  M=0.9  solution  by  the  procedure  of  Case  V is 
well  worth  the  extra  care  that  must  be  exercised  in  the  selection  of 
acceleration  parameters.  Prior  knowledge  of  the  effects  of  acceleration 
parameters  on  the  acquisition  of  an  acceptable  solution  is  essential. 

Case  VI  represents  an  Mot=0.9  condition  with  a strong  shock  in  the 
boattail  region.  The  pressure  distribution  over  the  body  is  illustrated 
in  Figure  26  . The  maximum  Mach  number  before  the  shock  is  1.295.  By 
using  an  optimum  AER  of  1.2  and  ABC=AHR=0.85  the  solution  converged 
smoothly  in  187  iterations.  Although  the  error  reached  the  convergence 
criteria  when  AER  was  set  to  1.25,  the  maximum  error  was  fluctuating 
as  it  did  so  on  the  202nd  iteration.  Changing  AHR  and  ABC  had  very 
little  effect  on  the  convergence  until  they  caused  divergence  when  set 
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to  1.0.  The  only  effect  shown  was  that  convergence  did  require  197 
iterations  when  AHR  was  reduced  to  0.60  with  AER=1 . 2 and  ABC=0.85.  A 
good  characteristic  that  this  solution  did  show,  and  which  the  other 
five  cases  lacked,  was  a satisfactory  flow  development  with  less  than  optimum 
values  of  AER. 

3.  BOATTAIL  PRESSURE  DRAG  CALCULATIONS 

A study  of  boattail  pressure  drag  calculated  by  the  inviscid  flow 
analysis  of  this  report  was  made  using  the  experimental  data  of  Reubush 
(Reference  17).  The  present  method  used  to  calculate  pressure  drag  is 
discussed  in  Section  III  D. 

Reubush  ran  several  cone-cylinder,  circular-arc  boattail  bodies 
with  solid  sting  plume  simulators  in  the  NASA  LaRC  16T  wind  tunnel. 

The  maximum  diameter  of  his  body  was  0.5  feet.  The  cone-cylinder  fore- 
body arrangement  was  4.0  feet  to  the  start  of  the  boattail.  The  flow 
region  was  simulated  with  a coordinate  system  similar  to  Figure  17. 

The  upper  boundary  was  set  at  the  proper  height  to  simulate  the  radius 
of  the  octagonal  16T  test  section.  A solid-wall  boundary  condition  was 
used  to  simulate  the  wall  even  though  16T  has  slotted  walls.  Two 
configurations  were  studied. 

The  first  (illustrated  in  Figure  25)  had  a boattail  length  to 
maximum  body  diameter  ratio  (L/D)  of  1.768.  This  boattail  had  a small 
separated  flow  region  near  the  junction  with  the  sting  which  starts 
at  about  station  1.6,  as  shown  in  Figure  25.  Comparisons  were  made 
at  flow  Mach  numbers  of  0.4,  0.6,  0.8,  0.9,  and  0.94.  The  present 
calculations  are  typical  of  inviscid  boattail  solutions  (Figure  25). 

In  the  region  ahead  of  the  separation  point,  the  pressure  distributions 
compare  very  well.  We  would  expect  the  drag  integrations  to  this  point 
to  compare  favorably.  However,  after  this  point,  the  inviscid  solution 
continues  to  recompress  while  the  actual  flow  separates  and  fails  to 
recompress.  The  result  is  a lower  drag  prediction  from  the  inviscid  analy- 
sis. This  fact  is  demonstrated  by  the  drag  plots  of  Figure  27.  All  across 
the  Mach  number  range,  the  inviscid  drag  is  lower  than  the  experimental 
drag. 
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X/0(t1AX)  - HEASUREO  FROM  START  OF  BOATTAIL  AT  BODY  STATION  X/0=8 


Figure  25.  Pressure  Distributions  for  Flow  over  the  NASA  Langley 
Circular  Arc  Boattail  (L/D  = 0.8)  at  M=0.4,  0.6,  0.8, 
0.9,  and  0.94 
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Figure  25.  (Continued) 
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X/OIMAX)  - MEASURED  FROM  START  OF  0OATTAIL  AT  BOOT  STATION  X/0=8 

Figure  25.  (Concluded) 
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X/0( MOX ) - HEASUREO  FROM  START  OF  BOATTAIL  AT  BOOY  STATION  X/0=8 


Figure  26.  Pressure  Distributions  for  Flow  over  the  NASA  Langley 
Circular  Arc  Boattail  (L/D  = 0.8)  at  M=0.4,  0.6,  0.8, 
and  0.9 
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Figure  27.  Pressure  Drag  Comparisons  for  NASA  Langley 
Circular  Arc  Boattails 
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The  second  configuration  featured  a much  shorter  boattail  with 
L/D=0.8.  It  was  expected  that  the  inviscid  analysis  would  not  predict 
the  pressure  distribution  over  this  highly  viscous  influenced  boattail 
flowfield  very  well.  Comparisons  were  made  at  Mach  numbers  of  0.4,  0.6, 
0.8,  and  0.9.  The  present  results  and  Reubush's  experimental  data  are 
plotted  in  Figure  26.  For  M=0.4,  the  pressure  distributions  compare 
well  up  to  the  separation  point.  From  there  on,  the  experimental  data 
levels  out  over  the  separation  region  and  the  inviscid  analysis  rises 
to  a typical  large  recompression  value.  As  expected,  the  drag  predic- 
tions by  the  invicid  analysis  are  lower  than  the  values  calculated  from 
the  experimental  pressure  distributions.  The  same  trend  continues 
through  M=0.8.  The  calculations  overshoot  the  experimental  pressures 
upstream  of  the  separation  point  and  cause  a calculated  drag  increase; 
the  calculated  recompression  error  grows  and  causes  an  even  larger 
calculated  drag  reduction.  At  M=0.9,  a shock  has  formed  on  the  boat- 
tail,  and  due  to  the  large  viscous  effects  the  calculation  does  not 
compare  well  aft  of  X/D(MAX)>0.3.  However,  note  that  drag  rise  is 
predicted  by  the  invicid  analysis  as  drag  increases  drastically  near 
M=0.9. 

' 

Although  the  invicid  drag  values  do  not  predict  the  exact  magnitude 
of  the  experimental  drag,  the  method  does  predict  the  drag  trends  and 
drag  rise  with  reasonable  accuracy. 

4.  UPPER  BOUNDARY  CONDITION  EFFECTS 

Experimental  studies  have  been  made  with  various  testing  techniques 
to  examine  the  performance  of  aircraft  nozzle  afterbodies  (NAB).  These 
testing  techniques  involved  the  application  of  various  upper  boundary 
conditions  to  the  fluid  flow  problem  as  well  as  different  support  and 
jet  plume  simulation  methods.  The  present  method  is  used  to  examine 
the  trends  exhibited  by  two  of  these  testing  techniques. 

First,  a description  of  the  two  test  configurations  is  in  order. 
Bittrick  (Reference  18)  made  a study  of  NAB  drag  using  a free  jet  model 
technique.  The  General  Dynamics  Free  Jet  Duct  Calibration  Facility  was 
used  with  the  NAB  attached  to  a strut- supported  hollow  shaft  that 
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located  the  model  centrally  in  the  9.25  inch  square  test  section.  Jet 

flow  to  the  nozzle  was  supplied  through  the  shaft.  The  boundary  layer 

entering  the  NAB  region  on  the  model  developed  along  a 4.0  foot  growth 

run  at  R =5.5  X 10®.  The  maximum  diameter  of  the  model  was  3.55  inches, 
e 

The  model  extended  9.05  inches  past  the  end  of  the  square  test  section. 

Bowers  (Reference  19)  tested  a similar  NAB  contour,  using  an  annular 
blowing  sting  technique.  The  28.68  inch  cone-cyl inder-boattail  model 
was  sting-mounted  in  the  Air  Force  Flight  Dynamics  Laboratory  Trisonic 
Gasdynamics  Facility.  The  facility  is  a closed  circuit,  variable  pres- 
sure wind  tunnel  with  a solid  wall  2.0  foot  square  test  section.  The 
jet  flow  was  supplied  through  the  hollow  sting  and  blown  back  annularly 
around  the  sting  to  simulate  the  nozzle  exhaust.  At  the  M=0.8  test 

condition  considered,  R =2.0  X 10®.  The  maximum  diameter  of  the  model 

e 

is  3.75  inches,  and  the  length  of  the  NAB  is  5.02  inches.  Bowers  ob- 
tained his  NAB  coordinates  from  a blueprint  of  Bittrick's  model. 

A study  using  the  analysis  described  in  Section  III  was  made  on 
three  upper  boundary  configurations  at  M=0.8.  The  coordinates  for  the 
calculation  were  obtained  from  construction  specification  dimensions 
from  Bowers'  NAB  test.  The  nominal  upper  boundary  condition  for  the 
potential  flow  is  a simulated  free  air  boundary  calculated  by  placing 
a solid-wall  boundary  condition.  Equation  49,  at  2.0  body  lengths 
(57.36  inches)  radius.  Even  with  this  free-air  boundary,  the  flow  over 
the  mid  section  of  the  body  did  not  quite  reach  the  free-stream  value. 

At  body  station  17.0  inches,  the  minimum  Mach  number  on  the  body  was 
calculated  to  be  0.806.  This  theoretical  result  indicates  that  the 
fineness  ratio  of  the  body  (7.6)  was  not  quite  large  enough  to  ensure 
that  the  approach  Mach  number  to  the  NAB  was  independent  of  forebody 
shape.  The  NAB  surface  pressure  distribution  for  the  free-air  upper 
boundary  case  is  plotted  in  Figure  28  and  flowfield  Mach  number  contours 
are  plotted  in  Figure  29. 

The  next  upper  boundary  condition  calculated  was  designed  to  simu- 
late the  AFFDL  2.0  foot  wind  tunnel  configuration.  Here,  as  with  the 
nominal  condition,  the  exhaust  plume  was  simulated  with  a cylinder  of 
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Figure  29.  Boundary  Condition  Effects  on  Nozzle  Afterbody 
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the  same  diameter  as  the  NAB  jet  exit.  Both  plots  are  again  presented 
in  Figures  28  and  29.  The  calculated  pressure  distribution  Indicates 
that  the  solid-wall  upper  boundary  at  1.0  foot  radius  from  the  center- 
line  forces  the  flow  to  accelerate  slightly  over  the  body.  The  minimum 
Mach  number  on  the  body  was  0.814  which  occurred  at  the  17.7  Inch  body 
station.  The  approach  Mach  number  to  the  NAB  was  correspondingly  higher. 

The  final  upper  boundary  condition  calculated  was  modeled  after  the 
free-jet  setup  of  Bittrlck.  This  calculation  was  made  with  no  forebody. 
The  upstream  boundary  condition  of  a uniform  stream  at  M=0.8  was  applied 
at  an  effective  body  station  of  10.0  inches.  The  free-jet  boundary  was 
approximated  by  placing  a straight  horizontal  boundary  at  R/D(MAX)=1 .305. 
The  free-jet  boundary  has  been  approximated  by  Murman  (Reference  20)  in 
small  perturbation  theory  by  setting  the  perturbation  velocity  equal  to 
zero.  For  this  method,  the  corresponding  treatment  was  to  solve  u^ 
along  the  free-jet  boundary.  The  NAB  coordinates  were  the  same  as  those 
used  In  the  other  calculations.  The  calculated  Mach  number  at  the  17.0 
Inch  body  station  Is  0.80001.  Therefore,  the  approach  Mach  number  to 
the  NAB  is  very  close  to  the  desired  value  of  0.8.  The  calculated  pres- 
sure distribution  and  Mach  number  contours  are  presented  in  Figures  28 
and  29. 

A comparison  of  the  calculated  flowflelds  may  be  made  by  consider- 
ing the  Mach  number  contour  plots  in  Figure  29.  Note  that  the  place- 
ment of  the  solid  wall  at  the  12.0  inch  position  causes  a growth  of 
the  supersonic  pocket  at  the  shoulder  and  the  appearance  of  a new  M=0.9 
contour  In  the  NAB  region  as  compared  to  the  free-air  case.  The  free- 
jet  case  predicts  a further  reduction  in  the  size  of  the  contours. 

In  addition  to  the  calculated  surface  pressure  distributions 
plotted  In  Figure  28,  the  experimental  data  of  Bittrlck  and  Bowers 
are  plotted  for  M=0.8  and  NPR=2.0.  It  is  felt  that  comparisons  of 
the  experimental  and  analytical  pressure  distributions  are  only  valid 
over  the  first  part  of  the  NAB  as  viscous  smoothing  effects  and  the 
imprecise  body  coordinate  definition  are  considerable  factors  in 
transonic  NAB  flow.  It  can  be  seen  that  the  calculated  pressures 
for  the  2.0  foot  configuration  and  Bowers'  data  compare  well  near  X/L=1.0. 
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Also,  the  level  of  the  calculated  pressures  for  the  free-jet  case 
compare  well  with  Bittrick's  data  for  X/L>0.6,  except  for  the  dip 
in  the  calculated  pressure  at  X/L=0.95.  It  is  expected  that  this  is 
due  to  a combination  of  viscous  smoothing  and  an  actual  difference 
in  NAB  coordinates  at  the  position. 

The  above  analysis  has  shown  that  the  present  analytical  method 
is  capable  of  correctly  handling  various  upper  boundary  conditions 
to  aid  in  the  interpretation  of  NAB  model  test  data  from  various 
techniques. 
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SECTION  V 
SUMMARY 

The  exact  velocity  potential  equation  for  steady,  two-dimensional 
and  axi symmetric  flow  in  simply  connected  regions  has  been  solved  with 
the  use  of  automatically  generated  body-fitted  curvilinear  coordinates 
from  incompressible  conditions  to  Mach  numoers  near  unity.  The 
following  is  a summary  of  the  results  of  this  effort: 

1.  The  definitions  of  the  attraction  factors,  P and  Q,  as  devel- 
oped by  Thompson  (Reference  1)  have  been  improved  to  ensure  predictable 
results  when  generating  coordinates  that  involve  attraction  to  both 
upper  and  lower  boundaries. 

2.  A new,  fast  system  of  generating  near  orthogonal  coordinates 
for  moderately  thick  (D/L < 0.3)  bodies  has  been  developed  and  proven 
to  be  reliable. 

3.  The  potential  solution  developed  has  been  shown  to  be  accurate 
and  reasonably  efficient  for  a wide  variety  of  body  shapes,  boundary 
conditions,  and  Mach  numbers. 

4.  An  axisymmetric  body  pressure  drag  integration  method  has  been 
applied  and  has  been  shown  capable  of  accurately  predicting  drag  rise 
and  drag  trend. 

5.  The  implementation  of  various  boundary  conditions  has  proved 
it  possible  to  simulate  various  wind-tunnel  testing  techniques. 
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